
exp1 <- read.table(file="exp_noed_nohet_randomcens0.txt", header=TRUE)
exp2 <- read.table(file="exp_noed_nohet_randomcens0.1.txt", header=TRUE)
exp3 <- read.table(file="exp_noed_nohet_randomcens0.2.txt", header=TRUE)
exp4 <- read.table(file="exp_noed_nohet_randomcens0.3.txt", header=TRUE)
exp5 <- read.table(file="exp_noed_nohet_randomcens0.4.txt", header=TRUE)
exp6 <- read.table(file="exp_noed_nohet_randomcens0.5.txt", header=TRUE)


betafrailgap <- array(NA, dim=c(6,1))
betafrailgap[1,] <- mean(exp1$betafrailgap, na.rm=TRUE)
betafrailgap[2,] <- mean(exp2$betafrailgap, na.rm=TRUE)
betafrailgap[3,] <- mean(exp3$betafrailgap, na.rm=TRUE)
betafrailgap[4,] <- mean(exp4$betafrailgap, na.rm=TRUE)
betafrailgap[5,] <- mean(exp5$betafrailgap, na.rm=TRUE)
betafrailgap[6,] <- mean(exp6$betafrailgap, na.rm=TRUE)

betafrailelapse <- array(NA, dim=c(6,1))
betafrailelapse[1,] <- mean(exp1$betafrailelapse, na.rm=TRUE)
betafrailelapse[2,] <- mean(exp2$betafrailelapse, na.rm=TRUE)
betafrailelapse[3,] <- mean(exp3$betafrailelapse, na.rm=TRUE)
betafrailelapse[4,] <- mean(exp4$betafrailelapse, na.rm=TRUE)
betafrailelapse[5,] <- mean(exp5$betafrailelapse, na.rm=TRUE)
betafrailelapse[6,] <- mean(exp6$betafrailelapse, na.rm=TRUE)

betaag <- array(NA, dim=c(6,1))
betaag[1,] <- mean(exp1$betaag, na.rm=TRUE)
betaag[2,] <- mean(exp2$betaag, na.rm=TRUE)
betaag[3,] <- mean(exp3$betaag, na.rm=TRUE)
betaag[4,] <- mean(exp4$betaag, na.rm=TRUE)
betaag[5,] <- mean(exp5$betaag, na.rm=TRUE)
betaag[6,] <- mean(exp6$betaag, na.rm=TRUE)

betacondgap <- array(NA, dim=c(6,1))
betacondgap[1,] <- mean(exp1$betacondgap, na.rm=TRUE)
betacondgap[2,] <- mean(exp2$betacondgap, na.rm=TRUE)
betacondgap[3,] <- mean(exp3$betacondgap, na.rm=TRUE)
betacondgap[4,] <- mean(exp4$betacondgap, na.rm=TRUE)
betacondgap[5,] <- mean(exp5$betacondgap, na.rm=TRUE)
betacondgap[6,] <- mean(exp6$betacondgap, na.rm=TRUE)

betacondelapse <- array(NA, dim=c(6,1))
betacondelapse[1,] <- mean(exp1$betacondelapse, na.rm=TRUE)
betacondelapse[2,] <- mean(exp2$betacondelapse, na.rm=TRUE)
betacondelapse[3,] <- mean(exp3$betacondelapse, na.rm=TRUE)
betacondelapse[4,] <- mean(exp4$betacondelapse, na.rm=TRUE)
betacondelapse[5,] <- mean(exp5$betacondelapse, na.rm=TRUE)
betacondelapse[6,] <- mean(exp6$betacondelapse, na.rm=TRUE)


censrate <- array(NA, dim=c(6,1))
censrate[1,] <- mean(exp1$censrate, na.rm=TRUE)
censrate[2,] <- mean(exp2$censrate, na.rm=TRUE)
censrate[3,] <- mean(exp3$censrate, na.rm=TRUE)
censrate[4,] <- mean(exp4$censrate, na.rm=TRUE)
censrate[5,] <- mean(exp5$censrate, na.rm=TRUE)
censrate[6,] <- mean(exp6$censrate, na.rm=TRUE)

yvalues <- c(betafrailgap, betacondgap, betacondelapse, betaag, betafrailelapse)
xvalues <- c(rep(censrate,5))
names <- c(rep("Conditional Frailty",6), rep("Conditional, Gap",6), rep("Conditional, Elapse",6), rep("Andersen-Gill",6), rep("Frailty, Elapse",6))
noednohet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

#library(ggplot2)
#betaplot1 <- ggplot(, aes(censrate)) + geom_line(aes(y=betafrailgap, colour="Conditional Frailty"))  + geom_line(aes(y=betacondgap, colour="Conditional, Gap"))  + geom_line(aes(y=betacondelapse, colour="Conditional, Elapsed"))  + geom_line(aes(y=betaag, colour="Andersen-Gill"))  + geom_line(aes(y=betafrailelapse, colour="Frailty, Elapse"))+ scale_colour_hue("Model") + labs(y="Mean Beta", x="Censoring Rate")

exp1 <- read.table(file="exp_ed_nohet_randomcens0.txt", header=TRUE)
exp2 <- read.table(file="exp_ed_nohet_randomcens0.1.txt", header=TRUE)
exp3 <- read.table(file="exp_ed_nohet_randomcens0.2.txt", header=TRUE)
exp4 <- read.table(file="exp_ed_nohet_randomcens0.3.txt", header=TRUE)
exp5 <- read.table(file="exp_ed_nohet_randomcens0.4.txt", header=TRUE)
exp6 <- read.table(file="exp_ed_nohet_randomcens0.5.txt", header=TRUE)

betafrailgap <- array(NA, dim=c(6,1))
betafrailgap[1,] <- mean(exp1$betafrailgap, na.rm=TRUE)
betafrailgap[2,] <- mean(exp2$betafrailgap, na.rm=TRUE)
betafrailgap[3,] <- mean(exp3$betafrailgap, na.rm=TRUE)
betafrailgap[4,] <- mean(exp4$betafrailgap, na.rm=TRUE)
betafrailgap[5,] <- mean(exp5$betafrailgap, na.rm=TRUE)
betafrailgap[6,] <- mean(exp6$betafrailgap, na.rm=TRUE)

betafrailelapse <- array(NA, dim=c(6,1))
betafrailelapse[1,] <- mean(exp1$betafrailelapse, na.rm=TRUE)
betafrailelapse[2,] <- mean(exp2$betafrailelapse, na.rm=TRUE)
betafrailelapse[3,] <- mean(exp3$betafrailelapse, na.rm=TRUE)
betafrailelapse[4,] <- mean(exp4$betafrailelapse, na.rm=TRUE)
betafrailelapse[5,] <- mean(exp5$betafrailelapse, na.rm=TRUE)
betafrailelapse[6,] <- mean(exp6$betafrailelapse, na.rm=TRUE)

betaag <- array(NA, dim=c(6,1))
betaag[1,] <- mean(exp1$betaag, na.rm=TRUE)
betaag[2,] <- mean(exp2$betaag, na.rm=TRUE)
betaag[3,] <- mean(exp3$betaag, na.rm=TRUE)
betaag[4,] <- mean(exp4$betaag, na.rm=TRUE)
betaag[5,] <- mean(exp5$betaag, na.rm=TRUE)
betaag[6,] <- mean(exp6$betaag, na.rm=TRUE)

betacondgap <- array(NA, dim=c(6,1))
betacondgap[1,] <- mean(exp1$betacondgap, na.rm=TRUE)
betacondgap[2,] <- mean(exp2$betacondgap, na.rm=TRUE)
betacondgap[3,] <- mean(exp3$betacondgap, na.rm=TRUE)
betacondgap[4,] <- mean(exp4$betacondgap, na.rm=TRUE)
betacondgap[5,] <- mean(exp5$betacondgap, na.rm=TRUE)
betacondgap[6,] <- mean(exp6$betacondgap, na.rm=TRUE)

betacondelapse <- array(NA, dim=c(6,1))
betacondelapse[1,] <- mean(exp1$betacondelapse, na.rm=TRUE)
betacondelapse[2,] <- mean(exp2$betacondelapse, na.rm=TRUE)
betacondelapse[3,] <- mean(exp3$betacondelapse, na.rm=TRUE)
betacondelapse[4,] <- mean(exp4$betacondelapse, na.rm=TRUE)
betacondelapse[5,] <- mean(exp5$betacondelapse, na.rm=TRUE)
betacondelapse[6,] <- mean(exp6$betacondelapse, na.rm=TRUE)


censrate <- array(NA, dim=c(6,1))
censrate[1,] <- mean(exp1$censrate, na.rm=TRUE)
censrate[2,] <- mean(exp2$censrate, na.rm=TRUE)
censrate[3,] <- mean(exp3$censrate, na.rm=TRUE)
censrate[4,] <- mean(exp4$censrate, na.rm=TRUE)
censrate[5,] <- mean(exp5$censrate, na.rm=TRUE)
censrate[6,] <- mean(exp6$censrate, na.rm=TRUE)

yvalues <- c(betafrailgap, betacondgap, betacondelapse, betaag, betafrailelapse)
xvalues <- c(rep(censrate,5))
names <- c(rep("Conditional Frailty",6), rep("Conditional, Gap",6), rep("Conditional, Elapse",6), rep("Andersen-Gill",6), rep("Frailty, Elapse",6))
ednohet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

#betaplot2 <- ggplot(, aes(censrate)) + geom_line(aes(y=betafrailgap, colour="Conditional Frailty"))  + geom_line(aes(y=betacondgap, colour="Conditional, Gap"))  + geom_line(aes(y=betacondelapse, colour="Conditional, Elapsed"))  + geom_line(aes(y=betaag, colour="Andersen-Gill"))  + geom_line(aes(y=betafrailelapse, colour="Frailty, Elapse"))+ scale_colour_hue("Model") + labs(y="Mean Beta", x="Censoring Rate")

exp1 <- read.table(file="exp_noed_het_randomcens0.txt", header=TRUE)
exp2 <- read.table(file="exp_noed_het_randomcens0.1.txt", header=TRUE)
exp3 <- read.table(file="exp_noed_het_randomcens0.2.txt", header=TRUE)
exp4 <- read.table(file="exp_noed_het_randomcens0.3.txt", header=TRUE)
exp5 <- read.table(file="exp_noed_het_randomcens0.4.txt", header=TRUE)
exp6 <- read.table(file="exp_noed_het_randomcens0.5.txt", header=TRUE)

betafrailgap <- array(NA, dim=c(6,1))
betafrailgap[1,] <- mean(exp1$betafrailgap, na.rm=TRUE)
betafrailgap[2,] <- mean(exp2$betafrailgap, na.rm=TRUE)
betafrailgap[3,] <- mean(exp3$betafrailgap, na.rm=TRUE)
betafrailgap[4,] <- mean(exp4$betafrailgap, na.rm=TRUE)
betafrailgap[5,] <- mean(exp5$betafrailgap, na.rm=TRUE)
betafrailgap[6,] <- mean(exp6$betafrailgap, na.rm=TRUE)

betafrailelapse <- array(NA, dim=c(6,1))
betafrailelapse[1,] <- mean(exp1$betafrailelapse, na.rm=TRUE)
betafrailelapse[2,] <- mean(exp2$betafrailelapse, na.rm=TRUE)
betafrailelapse[3,] <- mean(exp3$betafrailelapse, na.rm=TRUE)
betafrailelapse[4,] <- mean(exp4$betafrailelapse, na.rm=TRUE)
betafrailelapse[5,] <- mean(exp5$betafrailelapse, na.rm=TRUE)
betafrailelapse[6,] <- mean(exp6$betafrailelapse, na.rm=TRUE)

betaag <- array(NA, dim=c(6,1))
betaag[1,] <- mean(exp1$betaag, na.rm=TRUE)
betaag[2,] <- mean(exp2$betaag, na.rm=TRUE)
betaag[3,] <- mean(exp3$betaag, na.rm=TRUE)
betaag[4,] <- mean(exp4$betaag, na.rm=TRUE)
betaag[5,] <- mean(exp5$betaag, na.rm=TRUE)
betaag[6,] <- mean(exp6$betaag, na.rm=TRUE)

betacondgap <- array(NA, dim=c(6,1))
betacondgap[1,] <- mean(exp1$betacondgap, na.rm=TRUE)
betacondgap[2,] <- mean(exp2$betacondgap, na.rm=TRUE)
betacondgap[3,] <- mean(exp3$betacondgap, na.rm=TRUE)
betacondgap[4,] <- mean(exp4$betacondgap, na.rm=TRUE)
betacondgap[5,] <- mean(exp5$betacondgap, na.rm=TRUE)
betacondgap[6,] <- mean(exp6$betacondgap, na.rm=TRUE)

betacondelapse <- array(NA, dim=c(6,1))
betacondelapse[1,] <- mean(exp1$betacondelapse, na.rm=TRUE)
betacondelapse[2,] <- mean(exp2$betacondelapse, na.rm=TRUE)
betacondelapse[3,] <- mean(exp3$betacondelapse, na.rm=TRUE)
betacondelapse[4,] <- mean(exp4$betacondelapse, na.rm=TRUE)
betacondelapse[5,] <- mean(exp5$betacondelapse, na.rm=TRUE)
betacondelapse[6,] <- mean(exp6$betacondelapse, na.rm=TRUE)


censrate <- array(NA, dim=c(6,1))
censrate[1,] <- mean(exp1$censrate, na.rm=TRUE)
censrate[2,] <- mean(exp2$censrate, na.rm=TRUE)
censrate[3,] <- mean(exp3$censrate, na.rm=TRUE)
censrate[4,] <- mean(exp4$censrate, na.rm=TRUE)
censrate[5,] <- mean(exp5$censrate, na.rm=TRUE)
censrate[6,] <- mean(exp6$censrate, na.rm=TRUE)

yvalues <- c(betafrailgap, betacondgap, betacondelapse, betaag, betafrailelapse)
xvalues <- c(rep(censrate,5))
names <- c(rep("Conditional Frailty",6), rep("Conditional, Gap",6), rep("Conditional, Elapse",6), rep("Andersen-Gill",6), rep("Frailty, Elapse",6))
noedhet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)

#betaplot3 <- ggplot(, aes(censrate)) + geom_line(aes(y=betafrailgap, colour="Conditional Frailty"))  + geom_line(aes(y=betacondgap, colour="Conditional, Gap"))  + geom_line(aes(y=betacondelapse, colour="Conditional, Elapsed"))  + geom_line(aes(y=betaag, colour="Andersen-Gill"))  + geom_line(aes(y=betafrailelapse, colour="Frailty, Elapse"))+ scale_colour_hue("Model") + labs(y="Mean Beta", x="Censoring Rate")
#print(betaplot3)

exp1 <- read.table(file="exp_ed_het_randomcens0.txt", header=TRUE)
exp2 <- read.table(file="exp_ed_het_randomcens0.1.txt", header=TRUE)
exp3 <- read.table(file="exp_ed_het_randomcens0.2.txt", header=TRUE)
exp4 <- read.table(file="exp_ed_het_randomcens0.3.txt", header=TRUE)
exp5 <- read.table(file="exp_ed_het_randomcens0.4.txt", header=TRUE)
exp6 <- read.table(file="exp_ed_het_randomcens0.5.txt", header=TRUE)

betafrailgap <- array(NA, dim=c(6,1))
betafrailgap[1,] <- mean(exp1$betafrailgap, na.rm=TRUE)
betafrailgap[2,] <- mean(exp2$betafrailgap, na.rm=TRUE)
betafrailgap[3,] <- mean(exp3$betafrailgap, na.rm=TRUE)
betafrailgap[4,] <- mean(exp4$betafrailgap, na.rm=TRUE)
betafrailgap[5,] <- mean(exp5$betafrailgap, na.rm=TRUE)
betafrailgap[6,] <- mean(exp6$betafrailgap, na.rm=TRUE)

betafrailelapse <- array(NA, dim=c(6,1))
betafrailelapse[1,] <- mean(exp1$betafrailelapse, na.rm=TRUE)
betafrailelapse[2,] <- mean(exp2$betafrailelapse, na.rm=TRUE)
betafrailelapse[3,] <- mean(exp3$betafrailelapse, na.rm=TRUE)
betafrailelapse[4,] <- mean(exp4$betafrailelapse, na.rm=TRUE)
betafrailelapse[5,] <- mean(exp5$betafrailelapse, na.rm=TRUE)
betafrailelapse[6,] <- mean(exp6$betafrailelapse, na.rm=TRUE)

betaag <- array(NA, dim=c(6,1))
betaag[1,] <- mean(exp1$betaag, na.rm=TRUE)
betaag[2,] <- mean(exp2$betaag, na.rm=TRUE)
betaag[3,] <- mean(exp3$betaag, na.rm=TRUE)
betaag[4,] <- mean(exp4$betaag, na.rm=TRUE)
betaag[5,] <- mean(exp5$betaag, na.rm=TRUE)
betaag[6,] <- mean(exp6$betaag, na.rm=TRUE)

betacondgap <- array(NA, dim=c(6,1))
betacondgap[1,] <- mean(exp1$betacondgap, na.rm=TRUE)
betacondgap[2,] <- mean(exp2$betacondgap, na.rm=TRUE)
betacondgap[3,] <- mean(exp3$betacondgap, na.rm=TRUE)
betacondgap[4,] <- mean(exp4$betacondgap, na.rm=TRUE)
betacondgap[5,] <- mean(exp5$betacondgap, na.rm=TRUE)
betacondgap[6,] <- mean(exp6$betacondgap, na.rm=TRUE)

betacondelapse <- array(NA, dim=c(6,1))
betacondelapse[1,] <- mean(exp1$betacondelapse, na.rm=TRUE)
betacondelapse[2,] <- mean(exp2$betacondelapse, na.rm=TRUE)
betacondelapse[3,] <- mean(exp3$betacondelapse, na.rm=TRUE)
betacondelapse[4,] <- mean(exp4$betacondelapse, na.rm=TRUE)
betacondelapse[5,] <- mean(exp5$betacondelapse, na.rm=TRUE)
betacondelapse[6,] <- mean(exp6$betacondelapse, na.rm=TRUE)


censrate <- array(NA, dim=c(6,1))
censrate[1,] <- mean(exp1$censrate, na.rm=TRUE)
censrate[2,] <- mean(exp2$censrate, na.rm=TRUE)
censrate[3,] <- mean(exp3$censrate, na.rm=TRUE)
censrate[4,] <- mean(exp4$censrate, na.rm=TRUE)
censrate[5,] <- mean(exp5$censrate, na.rm=TRUE)
censrate[6,] <- mean(exp6$censrate, na.rm=TRUE)


#betaplot4 <-ggplot(, aes(censrate)) + geom_line(aes(y=betafrailgap, colour="Conditional Frailty"))  + geom_line(aes(y=betacondgap, colour="Conditional, Gap"))  + geom_line(aes(y=betacondelapse, colour="Conditional, Elapsed"))  + geom_line(aes(y=betaag, colour="Andersen-Gill"))  + geom_line(aes(y=betafrailelapse, colour="Frailty, Elapse"))+ scale_colour_hue("Model") + labs(y="Mean Beta", x="Censoring Rate")

yvalues <- c(betafrailgap, betacondgap, betacondelapse, betaag, betafrailelapse)
xvalues <- c(rep(censrate,5))
names <- c(rep("Conditional Frailty",6), rep("Conditional, Gap",6), rep("Conditional, Elapse",6), rep("Andersen-Gill",6), rep("Frailty, Elapse",6))
edhet <- data.frame(msqerror=yvalues,sample=xvalues,Model=names)


## B & W
library(ggplot2)

msqplot1 <- qplot(sample,msqerror,data=noednohet,geom="line",linetype=Model) + theme_bw()+  labs(y=expression(paste("Mean ",hat(beta),"  (",beta==-1,")")), x="Censoring Rate") + opts(legend.position=c(.7,.7)) + scale_linetype_manual("Model",values=c(1,2,4,3,5)) + opts(title="No Event Dependence, No Heterogeneity") +ylim(-1.3,-.7) 

msqplot2 <- qplot(sample,msqerror,data=ednohet,geom="line",linetype=Model) + theme_bw() +  labs(y=expression(paste("Mean ",hat(beta),"  (",beta==-1,")")), x="Censoring Rate")  + opts(legend.position="none") + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="Event Dependence, No Heterogeneity")  #+ylim(-.01,.6)

msqplot3 <- qplot(sample,msqerror,data=noedhet,geom="line",linetype=Model) + theme_bw() +   labs(y=expression(paste("Mean ",hat(beta),"  (",beta==-1,")")), x="Censoring Rate") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="No Event Dependence, Heterogeneity") 

msqplot4 <- qplot(sample,msqerror,data=edhet,geom="line",linetype=Model) + theme_bw()+  labs(y=expression(paste("Mean ",hat(beta),"  (",beta==-1,")")), x="Censoring Rate") + opts(legend.position="none")  + scale_linetype_manual("Model",values=c(1,2,4,3,5))  + opts(title="Event Dependence, Heterogeneity")


pdf("randomcens_panel.pdf",height=8.5,width=11)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2,2)))
vplayout <- function(x,y)
viewport(layout.pos.row=x,layout.pos.col=y)
print(msqplot1,vp=vplayout(1,1))
print(msqplot2,vp=vplayout(1,2))
print(msqplot3,vp=vplayout(2,1))
print(msqplot4,vp=vplayout(2,2))
dev.off()
                     
